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ABSTRACT 

A general formulation of the equivalent domain integral (EDI) method for 
mixed-mode fracture problems in cracked solids is presented. The method is dis- 
cussed in the context of a 3-D finite element analysis. The J-integral consists of 
two parts: the volume integral of the crack front potential over a torus enclosing 
the crack front and the crack surface integral due to the crack front potential plus 
the crack-face loading. In mixed-mode crack problems the total J-integral is split 
into and Jjjj representing the severity of the crack front in three modes 

of deformations. The direct and decomposition methods are used to separate the 
modes. These two methods were applied to several mixed-mode fracture problems 
in isotropic materials. Several pure and mixed-mode fracture problems were an- 
alyzed and results found to agree well with those available in the literature. The 
method lends itself to be used as a post-processing subroutine in a general purpose 
finite-element program. 


INTRODUCTION 

Several numerical techniques, in conjunction with finite-element (F-E) anal- 
yses, have been developed to calculate fracture mechanics parameters (stress- 
intensity factor K, strain energy release rate G , and J-integral). Three of these 
techniques are: (1) the virtual crack extension (VCE) method [1-4], (2) the virtual 
crack closure technique (VCCT) [5-8], and (3) the J-integral method [9-12]. The 
VCCT method is simple and accurate but can be applied only to linear elastic 
problems. In contrast, the VCE and J-integral methods can be applied to both 
linear and nonlinear problems. These methods are best demonstrated for pure 
mode problems or for calculating the total crack driving forces (G or J). Appli- 
cation of these methods to mixed-mode fracture problems is complex. The VCE 
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method involves a physical extension of the crack front by a small amount. Se- 
lection of the amount of crack extension is arbitrary and can introduce errors in 
inelastic problems. Further, both VCCT and VCE methods require a F-E mesh 
that is nearly normal to the crack front. Except for some simple cases, generating 
such customized F-E models for irregular-shaped cracks is difficult and time con- 
suming, if not impossible. Furthermore, such detailed modeling may not improve 
the global accuracy of the boundary value solution. Therefore, pursuit of methods 
that do not have these limitations continues. 

The J-integral method is very attractive, particularly for nonlinear material 
problems. With the original J-integral equation by Rice [9], Cherepanov [10,12], 
and Eshelby [11] for two-dimensional (2-D) problems as the starting point, several 
crack tip integrals were developed to include body forces due to thermal and 
magnetic fields and unloading effects in elastic-plastic problems [13-18]. For 2-D 
problems, the crack tip integrals are written as the sum of a remote line integral 
and an area integral around the crack tip. For 3-D problems, the J-integral is the 
sum of a remote surface integral and a volume integral around the crack front. 
These integral formulations suffer from a common drawback in that they require 
the evaluation of surface integrals which include singular terms. The evaluation 
of these surface integrals, although possible, is very unwieldy in F-E analyses. 

The J-integral formulation has been modified into a domain integral form [19- 
23] after de Lorenzi [24,25] introduced a S-function concept to define the virtual 
crack extension in 3-D cracked solids. In this method, the surface integrals for 3-D 
problems can be transformed into integrals over a domain or volume and, hence, 
the name equivalent domain integral (EDI). The EDI formulation is computation- 
ally very appealing and efficient. 

Recently, Nikishkov and Atluri [20] presented an EDI formulation for cracked 
3-D solids. To simulate the singularities at the crack front, they used quarter-point 
singularity elements. While they present the EDI formulation in a comprehensive 
manner, some details of the formulation need additional explanation. Also the 
formulation of reference 20 may not be general enough for problems where crack 
faces are subjected to external loading. The first objective of this paper is to 
present a general formulation of the EDI method for the calculation of J-integral 
under mixed-mode loading conditions. 
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Most investigators use collapsed quarter-point singularity elements at the 
crack front to simulate the crack front singularity with a polar arrangement of 
elements around the crack front. This type of mesh may be suitable but not con- 
venient, particularly, for crack extension studies. Furthermore, in the plastic range 
the quarter-point singular element produces a 1/r type singular strain field, which 
is valid only for elastic-perfectly-plastic material. Therefore, the second objective 
of this paper is to study the accuracy of the results when non-singular elements 
with a rectilinear arrangement of elements are used at the crack front. 

First the EDI formulation for general mixed-mode fracture problems in elastic, 
elastic-plastic, and anisotropic materials is presented. Next, the validity of the 
formulation is studied by applying it to several linear elastic and isotropic mixed- 
mode fracture problems. Numerical implementation of the EDI method for 20- 
noded and 8-noded, 3-D isoparametric elements is presented in the appendix. 
Several differences between the present formulation and those in the literature are 
highlighted. 


CRACK FRONT AND DOMAIN INTEGRALS 

The J-integral was introduced by Rice [9], Cherepanov [10], and Eshelby 
[11] to define the strength of the stress-strain field in nonlinear elastic 2-D crack 
problems. The J-integral was shown to be path-independent for nonlinear elastic 
and power-law hardening elastic-plastic materials. This path independence can 
be explained based on one singular point (crack tip) inside a closed contour in a 
singly connected domain. 

In 3-D crack problems, the crack front forms aline singularity and the strength 
of the singularity (K or J) could be varying all along the crack front. Therefore, the 
path independency is valid in a global sense, that is the total (or average) strength 
of the singularity of the complete crack front is independent of the surface enclosing 
it. However, at a point on the crack front, the path independency is valid only 
over a small region around the crack front due to interacting singular fields at 
neighboring points on the crack front. 

Consider a small tube of radius e around a segment of crack front of length 
A as shown in Figure 1(a) such that the limit of A and ^ tends to zero. The 
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local J-integral, also referred to as the crack front integral, over the surface A ( is 
defined as (see Fig. 1 (a) ) [12,20] 


J J Xh dx 3 = Um J[Wn k - dA ( 1 ) 

In Equation ( 1 ), W is the stress-work density, aij is the stress tensor, is the 
displacement vector, and Uk is the k th directional component of the unit normal 
vector on the closed surface A € . The indices i and j take the values 1,2, and 3, 
and k takes the values 1 and 2 . Thus the local value of J Xk is the total energy 
flux leaving the closed surface A e per unit crack front length in the k ih direction. 
J Xh can be defined in any coordinate system, but the local crack front coordinate 
system xi,Z 2 ? and x 3 is convenient for crack-extension studies. Note that the axes 
and x 3 are in the crack plane and are normal and tangential to the crack front, 
respectively, while x 2 is normal to the crack plane. 

The complete surface integral, in terms of surfaces identified in the Figure 
1 (a), is written as (henceforth the limits are dropped for convenience of presenta- 
tion) 


^ *7 Xh dx 3 


where 


f QdA + [ QdA + [ QdA 

Ja 4 J *fA <a J A« c| -M« cb 


■ 

Q = [Wn k - rij] 


( 2 ) 


W 


f€ij 

= I <7 ij deij 

Jo 


( 3 ) 

( 4 ) 


In Equation ( 2 ), A €l and A €2 are cross-sectional areas of the tube at 0 1 and O 2 , 
respectively. The subscripts ct and cb represent top and bottom crack surfaces, 
respectively. The stress-work density in Equation (4) is calculated over the com- 
plete strain path. The total strain tensor includes elastic, plastic, and thermal 
components. For linear elastic problems, W = - * Note that Equation ( 2 ) 

involves only stress, displacement and strain fields but no material properties. 
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However, to calculate stresses from strains the appropriate constitutive relation- 
ships (isotropic or anisotropic) must be used. Thus Equation (2) is applicable for 
general thermo-elastic-plastic and anisotropic material problems. 

Fracture modes in a 3-D cracked solid can be represented by three modes: 
opening (mode J), shearing (mode //), and tearing (mode III) modes. The 
corresponding three modes of the /-integral are J /, ///, and ////• In Equation 
( 1 ) or ( 2 ), J xl and J x2 represent the total /-integral (// + J// + Jm) and the 
product integral (-2y/Jj • ///), respectively [ 12 , 19, 20 ]. Since the meaning of J x3 
integral is not clear for crack problems, it is not defined by Equations (1) or (2). 
However, the mode III integral is separately defined as [20] 

f Jm dx 3 = lim [ f Qz dA + f Qz dA 

Jj\ ~k~*° J vA«i + A « 3 

A — ► 0 

+ f QzdA] (5) 

J A* ct + Atcb 

where 

Q 3 = W in ni - cr 3j rij (6) 

W 111 = f <Tij de 3 j (7) 

Jo 

The indices i and j take values 1, 2, and 3. Equations (2) through (7) completely 
define all three 7 -integrals. As previously mentioned, the evaluation of surface 
integrals, Equations (2) and (5), is tedious and could introduce errors due to 
numerical integration of singular terms. Therefore, an alternate form of evaluat- 
ing the above surface integrals called the Equivalent Domain Integrals (EDI) is 
presented in the next section. 

Equivalent Domain Integral 

Consider two tubular surfaces A ( and A spanning between two points O i and 
0 2 on the crack front (see Fig. 1 (b)). The tube A is arbitrary and encloses the 
tube A ( on which the 7 -integral is evaluated. The surface integrals (Eqs. ( 2 ) and 
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(5)) on A ( are converted into volume integrals by mathematical manipulations. 
First, the right hand side of Equations (2) and (5) are multiplied by unity. Next, 
from the resulting right hand side of the equation, subtract the product of the 
same integral over the surface A and zero, as shown below. 


L 


J Xk dx$ 


/ Jmdx 3 
Ja 


— 

1- 

f QdA 

v A f 

+ 1 • 

j QdA 

J A<i + A € 2 

+ 

1- 


QdA 

- 0 • J QdA 

= 

1. 

/ Q 3 dA 
Ja, 

+ 1 • 

f Q 3 dA 

•/ j4 € j -f A € 2 

+ 

1- 

*^A, ,, ■f A, c i> 

Qz dA 

— 0 * J Qz dA 


( 8 ) 


( 9 ) 


Note that, Equations ( 8 ) and (9) assume a unit extension of the crack front 
segment in the xj-direction. Instead, if an arbitrary (nonuniform) virtual extension 
of the crack front is made, Equations ( 8 ) and (9) need to be modified to account 
for this variation. Therefore, an arbitrary but continuous function S(x i, x 2 , x 3 ) is 
introduced [19, 20, 23] that has the property S(xi,x 2 ,X 3 ) = 0 on the surface A 
and S(xi,x 2 ,x 3 ) = S(x 3 ) on the surface A f . Using the S-function, Equations ( 8 ) 
and (9) are rewritten as follows. 


f J x ,Sdx 3 = [ QSdA + [ 

Ja J a € vA tl + a (J 


QSdA 


+ 


[ QSdA - f QSdA 

J A tci + a« c * Ja 


( 8 a) 


f JuiSdx 3 = f QzSdA + f Qz S dA 

Ja J At Ja € i + a € a 

Q 3 SdA - J Q 3 SdA 


+ 


L 


(9a) 


A 4c t + A € C 6 


As previously mentioned, A is small (lim A — * 0) and, hence, J Zh is assumed to 
be constant over the crack front segment length A . Then, Equations ( 8 a) and 
(9a) are simplified to 


6 



where 




f = f QSdA + f QSdA 

J A, J A, t + A, 2 


+ 


f QSdA - f QSdA 

J A €c t + A Mc b Ja 


Jiii • f 


= f Q 3 SdA + f Q 3 SdA 

J A 4 "At i + A €2 

+ [ Q 3 SdA - f Q 3 SdA 

J A €c * 4 - A* v A 


f 


B gt ^*cb 

»Oj 
'Ox 


r o 2 

= / S(x 3 ) 

JOx 


dx 3 


( 10 ) 

( 11 ) 

( 12 ) 


The parameter / is equivalent to the new crack surface area created by translating 
the crack front by 5 ( 213 ) in the xi- direction. Evaluation of / in Equation (12) 
and the choice of the 5-function will be discussed later. Equations (10) and 
(11) were further simplified by selecting the 5-function such that the function has 
zero values at two end surfaces {0\ and O 2 ) of the tubes A € and A and non-zero 
between these two end faces. With this choice, the second surface integrals in 
Equations (10) and (11) become identically zero. 

Now the integrals on the crack faces between the inner (A e ) and outer (-4) 
tubes are added and subtracted from the right hand sides of Equations (10) and 
(11). This manipulation is performed to obtain the integrals on a closed surface 
that encloses a volume. After some elementary algebraic operations, Equations 
(10) and (11) is rewritten as follows 


J Xk • / = - QSdA 

J A + (A — A,) el + A, + (A, — A) f t 

+ f QSdA + f QSdA (13) 

J(A — A,) ct J(A, — A) ch 

+ f QdA 

" A, c , + 
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Q 3 SdA 


Jm • / = - / 

J A 4- (A — -A«) cf + A t + (A € — A) cb 

+ f QzSdA + [ Q 3 SdA (14) 

J(A-A.) ct J(A.-A) cb 

+ f Q 3 dA 

J A, ct + -A« c » 

The first term in Equations (13) and (14) is negative because the direction of 
integration on the inner surface of the tube (.A*) is reversed. In Equations (13) 
and (14), ( A — A t ) ct and ( A ( — A) cb are the top and bottom crack surface areas 
between the two tubes A and A f . The first term in Equations (13) and (14) 
encloses the volume between the two tubes A and A e , which is represented as 
(V — V ( ). The rest of the terms in the right-hand sides of Equations (13) and 
(14) are integrals on the crack faces. Hence, J Xh and J/// are expressed as the 
sum of domain and crack surface integrals as follows 


J*>-f 

’ /) domain 

+ 

(^*1 1 * crack faces 

(15) 

Jill • f 

^ III * /) domain 

+ 

crack face* 

(16) 


Reference 20 obtained similar equations but did not include the crack face 
integrals. However, as will be shown later, for special cases the crack-face integrals 
vanish and the complete J-integral is given by the domain part of Equations (15) 
and (16). 

Domain integral .- Invoking Green’s divergence theorem, the closed surface 
integral of Equations (13) and (14) are written as a domain integral as follows 

in = ~j QSiA 


= -j 


\Wn k - SdA 


(17) 


J(v 


r d(WS ) d dui 
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Hence, the domain integral is 


(*^*» ’ f/ domain 


J(v 


r w ds dm dS 
[ W- — - (Tij- — -z— J dv 

— V € ) fc fc 


-l 


dW de {j 

[ ~ <Tij oT" ] 5 a*' 

(V-V.) 


Similarly one has 


[Jill * /) domain f t 

J(V- V € ) 




5a: i 


5ij 5xj 


-/ 


r 5W /JJ 5 ,5^ c JT/ 

[— ^i-5i-(«r:)] 5<iF 


(V-V.) 


dx\ 


dx\ v 5xj 


In deriving Equations (18) and (19), the equations of equilibrium 


( 18 ) 


(19) 


d(Ti 


*} 


dxj 


= 0 


and the small deformation strain-displacement relationships 


_ 1 , duj du j 

e,J 2 5x j dxi 


were used. 

In conventional finite element analysis, the equations of equilibrium are not 
satisfied point wise in the domain that is modeled. Numerical experimentation 
showed that the differences between including and not including the terms involv- 
ing the equations of equilibrium are of the order of 10 -3 to 10 -4 of the integral 
values for several problems. Therefore, in writing Equations (18) and (19) the 
equations of equilibrium are assumed to be satisfied exactly. 

The terms in brackets in the second integral in Equations (18) and (19) are 
point wise equal to zero for a linear elastic material. These terms, however, are 
non-zero in elastic-plastic and thermal problems. Since the present formulation is 
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for general thermo-elastic-plastic problems, such simplifications are not incorpo- 
rated. 

The domain integral Equation (18) is rewritten in a matrix form as 


(j.. ■/)<«*<. =- / [T#- {<J T k]{5’>]<iV 


dx k 

dW 


where 


- / i°} T U'*J)s dV 

J(V-V.) 9x k 

{<j} T = {(Til (722 033 0"l2 &23 031 } 


( 20 ) 


KJ T = 

r 5en 
1 5x* 

di-22 

dx k 

5633 

dxk 

9 5ei2 
5x* 

9 d £ 23 

dx k 

2 y ii} 

K> > T = 

5u : 
1 5x* 

0U2 

dx k 

du 3 

dx k 1 







<7X2 <731 






[*] = 

(712 

(722 <^23 







£731 

<723 <733 _ 






{5' } T = 

as 

{ 5xx 

as 

5X2 

55 

5x 3 ' 




(21) 


and the stress- work density W is 


W — I [<^ 11^11 4" (T22 d ^22 4“ 033^33 4“ 2<7]2(i£j2 4" 2^23^623 4" 2(73i<fe3i]. 

Jo 

Similarly, the mode III integral Equation (19) is written as 


— jUr n £- 


( 22 ) 


[ dWlH _ 

7(V-V.) ^*1 


^3} T {^'}]5 dF 
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where 


{* 3 } t 


{<T 3 i C 32 0 - 33 } 


f.nT _ r^3 d^U3 _&U 3_ 

1^3 / “loo O Q. o_ O- / 


and 


W 


III 


dx \ 2 0«i0*3 

= / [ 033^33 + *31 *31 + 0*32^32 ]• 

^0 


( 23 ) 


The numerical implementation of Equations (20) and (22) in a finite-element 
analysis with isoparametric elements is presented in the Appendix. 

The integral J Xl in Equation (20) for the linear elastic case is equivalent 
to the total strain-energy release rate calculated by the virtual crack extension 
method [1-4, 13, 24, 25]. 

Crack-face integrals .- The crack-face integrals in Equations (13) and (14) 

are 


(J Xh ■ f)crackface = [ Q S dA + f Q S dA 

— c t J(A € —A) ch 


^A 4< 


+ I QdA 

act A 4 c b 


(24) 


and 


( Jill ' /)cracfc/ac 


^(A — A t ) c 
^ A 4 c 


QsSdA + [ QzSdA 

J(A 4 — A)cb 


+ / Qz dA 

• ct cb 


(25) 


When the terms Q ( Eq. (3)) and Q 3 ( Eq. (6)) are zero on the crack faces, 
obviously, the integrals in Equations (24) and (25) vanish. On the crack faces, n 1 
and 713 are always zero, while 712 = —1 on the top face (ct) and 712 = 1 on the 
bottom face (cb). Imposition of these conditions in Equations (3) and (6), results 
in the following. 
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For k = 1 


(Q) 


crack face 


dui 

— <Ti2 -= — n 2 

OX I 


(,Q 3 )crackfa 


du 3 

— <T 32 q n 2 

OX i 


( 26 ) 


For k = 2 


(Q)cracJk/ace ^ ^^2 ^t2 ^ ^2 


(27) 


Note that for k = 2, Q3 is zero. 

Thus for traction free crack faces the terms (J Xl ) C rackface and (Jill) crack face 
vanish. In contrast, the term ( J X2 ) crack face is no longer zero. As noted in refer- 
ence 23, the ( J X3 ) C rackface is zero only for pure mode-7 fracture problems or for a 
singular stress field alone. However, in any finite size cracked body the stress field 
consists of both singular and non-singular terms and, hence, the (7 z 3 ) crack face in- 
tegral is not zero and cannot be neglected. The numerical evaluation of crack-face 
integrals involve the computation of singular integral terms, which are computa- 
tionally cumbersome and the source of errors [23]. 


Separation of Modes in Mixed-Mode Problems 

There are three modes of deformations in a cracked body, namely, the opening 
mode (Mode 7), the shearing mode (Mode 77), and the tearing mode (Mode 
777). The direct and decomposition methods are used to separate the mixed- 
mode fracture mechanics parameters into the three individual modes. 

Direct method .- The three components of the ./-integral, namely, 7/, 7// and 
7///, are calculated from J Xx and J Xj of Equation (15) and 7/jj of Equation 
(16). Since 7/// is directly calculated, the other two are calculated by solving the 
equations 
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Jl + Ju + Jill 


J X\ 

Jx 2 = —%V JlJlI ( 2 ®) 

Equation (28) was used to obtain J j and Jjj as 



Thus, in a general mixed-mode crack problem, computation of J Xl and J X7 
from Equation (13), Jui from Equation (14), and the use of Equation (29) 
completely defines all three modes of the J-integral. This procedure appears simple 
but the evaluation of J X7 could be complicated and erroneous due to the numerical 
integration of singular functions. Furthermore, as explained in reference 23, the 
local crack-face displacements are needed to distinguish between the opening and 
sliding (shearing) modes of deformation. Because of these reasons, separation of 
modes using Equation (29) may not be the best choice. Hence, an alternate 
method that avoids the evaluation of J X2 , called the decomposition method, is 
used. 

Decomposition method .- The advantage of transforming the surface inte- 
gral into a domain integral appears to be lost because of the non-zero crack-face 
integrals as shown in Equations (15), (16) and (27). These crack-face integrals are 
necessary to account for the terms containing the product of the singular and non- 
singular stress (strain) fields in the stress- work density expression. It was shown 
in reference 23 that the product terms are eliminated by decomposing the stress 
and displacement fields into symmetric and antisymmetric parts. The resulting 
equation contains only the domain integrals. Hence, the method is attractive 
and is computationally efficient. The decomposition method, however, requires 
additional effort to create a symmetric mesh about the crack plane. 
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The decomposition of displacement and stress fields is straight forward for 2- 
D problems [19, 23, 26-28], but is slightly more complicated for 3-D problems [20]. 
Hence, the decomposition of displacement and stress fields corresponding to the 
three modes of fracture is presented. Consider any two points P(x ux 2 ,xz) and 
P , (x 1 , — x 2 , xs ) that are in the immediate neighborhood of the crack front and are 
symmetric about the crack plane as shown in Figure 2. For any arbitrarily general 
deformation, the displacements and stresses at points P and P 1 can be expressed 
as a combination of symmetric and antisymmetric components as shown below. 


and 


' Uip ' 

f Uis' 



U 2 P 

> = < «2S 

> + < 

«2 AS 

< U *P J 

l U *S J 


[ «3 AS j 


{ U 1P > y 


UlS ) 


f -UlAS 

/ 

U2P> 

► = i 

~U2S 

► + < 

V-2AS 

[ u 3 p> J 

! 

{ U3S J 


{ -USAS J 


(30) 


(31) 


where subscripts S and AS denote the symmetric and antisymmetric components, 
respectively. 

Equations (30) and (31) are used to determine the symmetric and antisym- 
metric displacements in terms of the displacements at points P and P r as 


It 2 
^3 



tii 

ti 2 

uz 



1 

2 

1 

2 


. \ 
tiip + Uip* 

U2P — li2 P* > 
Uzp + tl3P' ; 

U\p — tiiP' 

U2P + ti2 P' > 

uzp — uzp> J 


(32) 


Similarly, the symmetric and antisymmetric components of the stresses are 
expressed in terms of the stresses at points P and P 1 (see Fig. 3) as 
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f o-n 1 


' <7ll p + <7ll p , ' 

<722 


<722f + <722,,, 

J <733 

1 

= - ( 

<733p + <733,. 

<712 

2 

<712,. — <712 p , 

<^23 


&23 P — <723pi 

k <731 > 

s 

k <73l p + <731^, > 


' <711 ' 


' <711 P ~ <7 lip, ’ 

<722 


<722 P ~ <722 p» 

C 33 

1 

= - ( 

<733 p ~ <733pr > 

<7i2 

2 

<7l2p + <7l2p f 

<723 


<723 p + <723 p r 

k <731 - 

AS 

k <731 p — <731 pi ' 


The symmetric and antisymmetric displacement and stress fields are further 
separated into mode J, mode II , and mode III components as follows 

M = {u} / + {u } 11 + {u } 111 


1 


U\P + Uipi 


= - < «2P - “2 P‘ 

u 3 p + U 3 pi 


1 

+ 2 < 


UiP — U\P> 
U 2 P + U 2 P> 
0 


1 

+ 2 


0 

0 

U 3 p — U 3 p> 


( 34 ) 


and 
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{»} = W' + M" + {«■}"' 

+ flip, 

cr 2 2 P + <7’22 f , 

033P + ViZpi 

<Tl2 P ~ <T 12p, 

CT23P — V2Z P , 

&3\p + Vzip, 

V\\p — CUpi 
<T 2 2p — 022p» 

0 

tr 12p + < 7 'l2p, 

0 
0 

0 
0 

^33 p — <733p, 

0 

^23p + ^23 p, 

^lp — ^lp, 

Similar equations in reference 20 had typographical errors. The mode I, 
mode II, and mode III displacements (Eq. (34)) and stresses (Eq. (35)) are 
used to directly evaluate Jj,Jij , and Jm from J Xl using Equation (15). Note 
that the surface integral in Equation (15) is required only when the crack face 
is loaded. The J Xi integral for each of these modes of deformation is identically 
equal to zero because of the orthogonality of the modes of deformations. Hence, 
the decomposition method involves only the evaluation of domain integrals and is 
computationally efficient. 

5-Functions 

As mentioned previously, the S-function is an arbitrary but continuous func- 
tion with a zero value on the surface A and at the ends of the tube ( A\ and A 2 ) 
and a non-zero value (varying between zero and one) on the surface A t (see Fig. 
1). On the tube surface A ( , the S-function is a function of only X 3 and has a value 
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of one at the location where the /-integral is required. The S-functions are con- 
veniently defined by specifying the values of S at the nodes and using the element 
shape functions (see Appendix). Figure 4 presents several types of S-functions for 
domains spanning one or two elements in the 13-direction. Table 1 presents the 
values of / (see Eq. (12)) that correspond to each of the S-functions shown in 
Figure 4. Also, the S-functions for both 8-noded and 20-noded isoparametric ele- 
ments are presented. For the 8-noded element, a linear S-function is defined, while 
for 20-noded elements several combinations of linear and quadratic functions are 
defined. Note that the values for / depend only on the variation of the S-function 
in the x 3 -direction. 

The /-integrals for various cracked 3-D solids were calculated using the various 
S-functions presented in Figure 4. The results of these numerical experiments will 
be discussed later. 


RESULTS AND DISCUSSION 

The EDI method was applied to several pure mode-/, mode-//, mode-///, and 
mixed-mode fracture problems to evaluate the accuracy, domain independency, 
and S-function independency. Although the method formulated above is for gen- 
eral anisotropic and nonlinear materials, the results presented here are restricted 
to linear elastic and isotropic materials with Poisson’s ratio of 0.3. Throughout the 
analysis only non-singular elements were used around the crack front. Wherever 
possible, a rectangular arrangement of finite elements was used near and around 
the crack front to evaluate the accuracy of the rectangular type of modeling. 

First, the EDI method was applied to a 3-D cracked body subjected to a com- 
bination of known mode-/, mode-//, and mode-/// singular stress fields. Then, 
the method was applied to various finite size crack problems subjected to loadings 
that produce single or mixed-mode deformations. The specimen configurations 
considered were the middle-crack tension specimen and embedded cracks in cir- 
cular cylindrical rods subjected to tension, torsion, and shearing loads. The com- 
puted /-integral values are compared with those from literature wherever possible. 
Both the direct and decomposition methods were used to separate the mode-/// 
component. Only the decomposition method was used to separate mode-/ and 
mode-// components because of the singular integration involved in the direct 
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method [23]. All of the above analyses used finite element models with 20-node 
isoparametric elements. Typical results are also presented for models with 8-node 
isoparametric elements. 

Singular Field Loading on a Cracked Body 

Consider a single-edge cracked specimen with a straight crack front as shown 
in the Figure 5(a). This crack configuration was subjected to mode-J, mode-// 
or mode-/// or a combination of these modes. This problem demonstrated the 
accuracy of the method without introducing the nonsingular stress field that occurs 
in any finite-element analysis of a cracked body. The displacements in each of the 
three modes of deformation are given in terms of the spherical coordinates (r, 0, 
andz) [29] as 
Mode-/ displacements: 



u 3 = 0 


Mode-// displacements: 



u 3 = 0 


Mode-/// displacements: 


(36) 


(37) 
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( 38 ) 


u\ = 0 
U2 = 0 


U3 



where G = E/{ 2(1 +*/)) 

The above displacement fields were imposed on the cracked body shown in Figure 
5(a). The resulting stresses and displacements in the solid were used to calculate 
the J-integrals by the EDI method. 

Utilizing the symmetries in the problem, a representative quarter of the solid 
was modeled with 20-node isoparametric elements as shown in Figure 5(b). One 
layer of elements was used to model the entire thickness. The finite element model 
had 320 nodes and 36 elements. The displacements for each mode of deformations 
were calculated at each node of the model from Equations (36)-(38) and were 
used as input for EDI algorithm. 

Domains and 5-functions .- Five domains were used in the calculations. 
These domains are shown in Figure 5(b) as D i , D 2 , . . ■ Ds • Each domain consisted 
of one ring of four elements around the crack front. The surface that is nearest to 
the crack front of each domain corresponds to the surface A f . Similarly, the surface 
that is farthest from the crack front of each domain corresponds to the surface A. 
(Thus, for example, the surface A for the domain D 2 will be the surface A ( for 
domain D 3 .) On the surfaces A of each domain, the 5-function was prescribed to 
be zero. On the surfaces A f several types of 5-functions in the X3- direction (the 
six types are shown in Fig. 4) were considered. Since one layer of elements was 
used to model the thickness of the solid, A 2 was set equal to zero for the Type II 
through VI 5-functions. 

■7-integral results .- Table 2 presents the normalized value of J Xl calculated 
for four domains (D 2 to D s ) and six types of 5-functions. The imposed displace- 
ment field on the model corresponds to K j = 1. The Jjj integral, as expected, 
was computed to the order of machine zero and hence is not shown here. Results 
for all four domains and Type J, II, III, and V 5-functions are in excellent 
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agreement with the exact solution. Type IV and VI 5-functions are quadratic in 
radial direction. Although these results are within two percent of the exact solu- 
tion they are not as accurate as other types. Similar trends were observed when 
quadratic 5-functions were used on several other crack problems. These results 
suggest that the simple Type-I 5-function will give accurate J-integrals. 

The computed integrals were inaccurate when the domain Di was used. The 
J Xl -integrals for this domain were about 10 to 15 percent higher than those for 
other domains. This behavior is attributed to the errors in the numerical integra- 
tion of the imposed singular stress field. 

Note that in the above analysis one ring of elements around the crack front 
wets used in each domain. Multiple rings of elements in each domain (for example, 
combining the Z> 2 and D 3 domains) gave results identical to those with a single 
ring of elements. Multiple rings, however, increase the data preparation efforts 
considerably. Thus, domains with a minimum number of rings of elements (in this 
case one ring) are preferable. 

The domain independency observed in this example is expected since this is 
a plane-strain problem where the strength of the singularity is constant along the 
crack front (i3-axis). Therefore, for crack bodies having a constant singularity 
strength along the crack front, the J Xl -integral is independent of the domain. 

Table 3 presents the normalized J Xl , calculated assuming several linear combi- 
nations of the Ki, K 11 , and Km displacement fields given by Equations (36-38). 
Four domains, each with Type I 5-functions were used. The J Xl calculated for 
four domains from Equation (15) agrees very well with the exact solutions. The 
maximum difference is less than 2% for the K 1 — Ku = 1 loading in domain D 3. 

Middle-Crack Tension Specimen 

A typical middle-crack tension, M(T)*, specimen of W/a = 2, t/a = 1, 
H / a = 8, and crack length a subjected to a uniform tension stress <r at 12 = i.H / 2 
was analyzed. The F-E mesh shown in the Figure 6 was used for the analysis by 
imposing symmetry conditions at X\ = —a, X2 = 0 (on the uncracked plane), and 
x 3 = 0 planes. The model was comprised of six unequal layers with thicknesses 
0.22t, 0.13t, 0.08t, 0.04t, 0.02t, and O.Olt as shown in Figure 6(a). (The layer with 

* ASTM abbreviation for middle-crack tension 
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the smallest thickness is near the 13 = t/2 surface.) Five domains, D 1, D 2 , • • • D 5, 
as shown in Figure 6(b) and Type I 5-function (linear in radial direction and 
quadratic in X3-direction, see Table 1) were used to calculate the domain integrals. 

Plane-strain analysis .- First, the M(T) specimen was analyzed for a plane- 
strain condition by imposing 113 = 0 on ij = 0 and 13 = t/2 planes. The J Xx 
values were calculated for the five domains all along the crack front using each 
of the six layers in the thickness direction. The average normalized value of J Xx , 
J Xl E/[ircr 2 a(l — i/ 2 )], was 1.410 with a maximum variation of less than 0.1%. 
This value agrees very well with handbook value of 1.414 and VCCT method [6] 
value of 1.424. These results also suggest that a simple linear 5-function in the 
radial direction yields accurate J Xl values. 

Three-dimensional analysis .- The M(T ) specimen was reanalyzed by re- 
laxing the plane-strain condition, i.e., by removing the boundary condition u 3 = 0 
on the X3 = t /2 plane. Three domain definitions, Da, Db, and Dc , were used. 
The domain and 5-function definitions for domain D 3 axe illustrated in the Figure 
6(b). In Da, each ring of elements around the crack front represents a domain. 
In this case, the radius (e) of the inner surface (the A e surface on which the J Xx 
is evaluated) was different for each domain; it varied from 0.0 to 0.734 times the 
crack length a for domains D\ to D 5 , respectively. The other two domains (Db 
and Dc) involve a constant inner surface A f and a variable outer surface. Domain 
Db had e = a/10, and, hence, the domain D 3 included the second and third rings 
of elements around the crack front. The corresponding 5-function definition in 
the radial direction is bilinear as shown in the Figure 6(b). The domain D c had 
c = 0. Therefore, for example, domain D 3 included first, second, and third rings 
of elements. 

Table 4 presents the normalized J Xl for all five domains using the three domain 
definitions. As expected, J Xl for the interior layer is independent of the domain 
and the domain definition. In contrast, J Xl for the last layer with a Da domain 
definition shows a strong domain dependency. As the radius (or mean distance 
from the crack front) of the inner surface of the domain becomes large the stress, 
strain, and displacement fields on this surface also include the effect of the singular 
field from the other segments of the crack front. Evaluation of J Xl -integral on this 
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surface would yield values which are different if the single crack front segment 
alone had contributed to the deformation field. If the inner surface of the domain 
is close to or at the crack front, the interaction from the neighboring segments of 
the crack is either very small or nonexistent. Hence, the Db (e/a = -1) an d Dc 
domain definitions gave more accurate results than the type Da- These results 
suggest that the radius of the inner surface should be less than or equal to one- 
tenth of the crack front. Furthermore, a domain consisting of only one ring of 
elements is sufficient to calculate accurate values of J provided the domain is very 
close to the crack front. Note that domains Di and D 2 always satisfy the above 
conditions, hence, the results were accurate and agreed very well with the VCCT 
method [6]. 

The global average value of J Xl over the crack front, however, is domain 
independent. Any ring of elements (even with a Da domain definition) spanning 
the complete length of the crack front yields accurate values. For example, the 
global average value of J Xl , E J Xl /[tt <t 2 a (1 — t' 2 )], for the four domains were 1.521 
( D 2 ), 1.523 ( D 3 ), 1.529 (D 4 ), and 1.531 (-D 5 )- Only the innermost domain was 
less accurate, 1.433 (Di ), because of errors in the stress-strain fields very near the 
crack front. In the above analyses a Da domain definition was used and each ring 
had 24 elements. 

Figure 7 compares the normalized J Xl along the crack front from EDI and 
VCCT [6] methods. Excellent agreement is observed between the two methods. 
Plane-strain results are also shown in the figure as a reference solution. 

Embedded Penny Shaped Crack in Circular Rod 

An embedded penny shaped crack of radius a in a circular cylindrical rod with 
D/a = 10 and H /a = 40 is shown in Figure 8. Two types of loadings, uniform 
tension and torsion, were considered. Note that x 3 , x 2 , and £3 represent the global 
coordinate system and u 2 , u 2 , and u 3 represent the corresponding displacements. 
Utilizing the symmetry in the problem one-eighth of the specimen was modeled. 
Figure 8(b) shows the F-E model of the lower eighth of the specimen. The model 
has 5143 nodes and 1020 twenty-noded elements with 6 layers in the circumferential 
direction as shown in Figure 8(b). Figure 8(c) shows the details of the modeling 
near a point on the crack front and the two domains D\ and D 2 used in EDI 
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calculations. Note the rectangular arrangement of elements. This arrangement 
is in contrast to the polar arrangement used in previous examples and by other 
investigators [20]. Two domains and Type I 5-functions were used to evaluate the 
integrals. 

Remote tension load .- The solid is loaded by a uniform stress <r on the 
x 3 = ±H /2 planes. The J Xl -integrals were calculated for the two domains (Da 
domain definition) in each of the six layers in the circumferential direction. The 
J xl values for all six layers calculated in either D\ or D 2 were identical. This 
is expected because of the axisymmetric nature of the problem. The normalized 
E J Xx l[cr 2 a(\ - i/ 2 )] values from the EDI method for domains D 1 and D 2 are 
1.333 and 1.299, respectively. The exact analytical solution due to Sneddon [30] 
for a penny shaped crack in an infinite solid is 1.273, while Benthem and Koiter s 
asymptotic solution [31] for a finite size solid is 1.275. Thus, the normalized J x \ 
values from domains D\ and D 2 differed by less than 3% from each other and 
are about 2% to 5% higher than the closed form solutions [30, 31]. The result 
for domain D 2 is more accurate than that for domain D\, about 2% higher than 
Benthem and Koiter’s value. The larger error in the domain D 1 may be attributed 
to the fact that only two elements were used around the crack front and the whole 
domain is within the singular field. A Dc domain definition of domain D 2 , which 
included the two rings of elements around the crack, also gave results identical to 
the D 2 domain solution. 

Remote torsion loading .- A torque of magnitude T = 7 t a*G/H was im- 
posed at the two ends of the specimen (13 = ±H/2), where G = E/[ 2(1 + v)}. 
This corresponds to an angular twist of magnitude 2 / H in the uncracked rod. The 
displacement field for this loading in the uncracked rod is 

til = —(2/H) x 2 

u 2 = ( 2 /IT)x 3 *i (39) 

U 3 = 0 

The boundary conditions u\ = 0 on x 2 =0 plane, u 2 = 0 on S\ =0 plane, 
and U 3 = 0 at all nodes were imposed on the finite-element model. On the face 
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X 3 = -H/2 the displacements given by Equation (39) were imposed on the 
model. The J — integral was calculated using both the direct and decomposition 
methods. The normalized values, E J jjj a 5 /[(l+i/) T 2 ], from the direct method for 
domains, £>i and D 2 are 0.2402 and 0.2389, and from the decomposition method 
for domains Z>j and D 2 are 0.2388 and 0.2374, respectively. The value of the 
normalized integral for an infinite solid obtained with the analytical solution by 
Lowengrub and Sneddon [32] was 0.2293. Since only mode - III loading was applied 
the total integral J Xl (Equation (15)) and the J/jj integral (Equation (16)) are 
nearly identical for both domains. The value of J jji calculated by the EDI method 
is about 4% larger than Lowengrub and Sneddon’s infinite body solution [32]. 


Inclined Embedded Penny Shaped Crack 


The EDI algorithm is next applied to problems involving mixed-mode defor- 
mations. Figure 9 shows an embedded inclined penny shaped crack in a circular 
rod subjected to a uniform stress a. The crack plane is inclined at an angle a to 
the ij — x 2 plane. The solution to this problem can be obtained as the sum of 
two solutions: the solution to a penny shaped crack subjected to traction normal 
to the crack faces of the magnitude cr n = (a/ 2)(1 + cos 2 a) and the solution 
of a penny shaped crack subjected to shear traction on the crack faces of magni- 
tude r = — (er/2) sin 2 a. The exact solutions for an infinite solid with the above 
mentioned traction are given by Cherepanov [12] for v = 0 and Kassir and Sih 
[33] for a general value of Poisson’s ratio v. The strain energy release rates for 
the three modes, after converting the stress-intensity factors given in reference 33 
using plane-strain assumptions, are 


Gi = 
Gu = 
Gin — 


a 2 a 
7 r E 


(1 + cos 2 a) 2 


4(1 — u 2 ) a a 2 . 2 2± 

7—— sm 2a cos a> 

7r(2 — v) 2 E 

4(1 — v)(l — u 2 )aa 2 

7T (2 — v) 2 E 


sin 2 2a sin 2 <f> 


(40) 


Where <f> is the angle measured from the zY-axis on the crack plane (see Fig. 9(b)). 
These results for the infinite size solid are used to compare with the results from 
EDI method for a finite size rod. 
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A cylindrical rod with aD/a = 5 and H/a = 10 was used in the analysis. Uti- 
lizing the symmetry in the problem one-half of solid was modeled. The model has 
9547 nodes and 2000 twenty-noded parabolic elements (see Fig. 10). Symmetry 
conditions (if 2 = 0) were imposed on x ‘2 = 0 plane. The three domains D 1 , D 2 , 
and Dz used are shown in Figure 10(b). The normalized total J-integral all along 
the crack front calculated using Equation (15) (J x 1 ) and the total value of J 
(sum of Jj, J/j, and Jjjj) from the decomposition method are shown in Figure 
11. The Kassir and Sih [33] solution is shown by the solid line in the figure. The 
total values of J from the direct and the decomposition methods agree with each 
other for all three domains. EDI results from all three domains are about 2% to 
7% larger than Kassir and Sih’s solution because a finite size solid was analyzed. 
Among the three domains, the domain D 2 g&ve the lowest value while the domain 
Z >3 gave the largest. The maximum difference between the two domains Z? 2 and 
Dz was about 3 percent. 

Figure 12(a) shows the distribution of the normalized Jj, J//, and Jm val- 
ues along the crack front obtained from the decomposition method. The indi- 
vidual mode components from the EDI method agree well with those determined 
by Kassir and Sih all along the crack front. The small differences between the 
two solutions can be attributed to the finite solid. It is interesting to note that 
Cherepanov’s solutions (with u = 0) for J jj and Jjji are about 28% lower than 
Kassir and Sih solution (with u = 0.3). 

Figure 12(b) compares the J m values calculated from the direct (Equation 
(16)) and decomposition methods with Kassir and Sih solution [33]. For all three 
domains, the two EDI procedures agree well with each other and with Kassir and 

Sih. 


Inclined Semi-Circular Surface Crack in a Tension Rod 


A semi-circular prismatic rod with an inclined semi-circular surface crack 
was analyzed. The centrally located semi-circular surface crack plane is oriented 
at an angle 30° to the xj-axis of the rod (see the insert sketch in Fig. 13). 
The rod is subjected to a remote tension loading. The finite-element model of 
Figure 10(a) was utilized to analyze this problem. The three domains D 1 , I? 2 , 
and Dz shown in Figure 10(b) were used in the EDI calculation. The total value 
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of J(J Xl ) and the individual components //, ///, and Jm along the crack front 
are shown in Figures 13 and 14, respectively. These results were obtained using 
the decomposition method. No reference solutions are available for this problem. 
Solid lines in Figures 13 and 14 are the smooth fit to the domain D 2 solutions. 

The /-values for the three domains agree well with one another. The maxi- 
mum difference between the results of the domain D 2 and D 3 is about 3 percent. 
As expected, (a) all the / values are symmetric about <f> = 90°, (b) the 7/ and J// 
have the highest values at </> = 0° and 180° with the lowest value at <f> = 90°, (c) 
the Jm value is zero at <f> = 0° and 180°, (d) Jm value is maximum at <f> = 90° 
and (e) J// value is zero at </> = 90°. 

Application to 8-Node Isoparametric Elements 

The use of 8-noded isoparametric elements in the EDI analysis was demon- 
strated using two classical crack problems in a tension specimen: an embedded 
elliptic crack and a penny shaped surface crack. The F-E model and the domain 
used were same as that given in reference 34, except the notch radius is zero. The 
F-E mesh had 3420 nodes and 2772 eight- noded elements. Around the crack 
front, a rectangular arrangement of non-singular elements was used (similar to 
that shown in Fig. 8(c)). The element size around the crack front was a/20, 
where a is the crack length in thickness direction of specimen. A domain consist- 
ing of the second ring of elements around the crack front and two element layers 
along the crack front was used for the /-integral calculation. 

Figure 15 compares the distribution of the normalized total J(J Xl ) along the 
crack front for an embedded elliptic crack calculated from the EDI method to 
Green and Sneddon’s infinite body solution [35]. The two solutions agree well 
with each other. 

Figure 16 shows the distribution of the normalized total / (which is same as 
7j) along the crack front for a semi-circular surface crack in a tension specimen. 
The /-values from the EDI method agree very well with those from the VCCT [6] 
and Raju and Newman’s solutions [36]. The EDI algorithm is incorporated in the 
ZIP3D [37] code. The ZIP3D is an elastic and elastic-plastic finite-element code 
to analyze cracked bodies. 
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CONCLUDING REMARKS 


Details of development of the three dimensional equivalent domain integral 
(3-D EDI) method for the calculation of mixed-mode fracture mechanics param- 
eters in isotropic or anisotropic and linear or nonlinear materials are presented. 
Differences and improvements between the current algorithm and that reported 
in the literature are highlighted. Results presented in this paper are restricted to 
isotropic elastic solutions. Several single and mixed-mode loaded cracked bodies 
were analyzed and results were found to agree very well with those available in 
the literature. 

The EDI method with 20- or 8-noded isoparametric, nonsingular elements 
and either a polar or rectilinear arrangement of elements at the crack front gives 
accurate values of the /-integral. A simple linear S-function in the radial direction 
is recommended if only one ring of elements at the crack front is used. The EDI 
method was found to be independent of the type of S-function, except for one 
special case. The / values were found to be inaccurate for domains consisting 
of only one ring of elements around the crack front with S-functions that were 
quadratic in the radial direction from the crack front. 

The EDI method is domain independent provided the radius of the inner 
surface of the domain is either zero or very small (less than one-tenth of the major 
crack length). Domains consisting of only the second ring of elements or the first 
two rings of elements around the crack front reduce the data preparation effort 
and also give accurate /-integral values. 

The principal advantage of the 3-D EDI method is that the finite element 
idealization need not be orthogonal to the crack front. The orthogonality of the 
modeling at the crack front is a requirement for the virtual crack closure and the 
force methods. In the case of mixed-mode loadings, the decomposition method 
yielded accurate /-integrals. The method requires the evaluation of only one 
integral with different sets of displacement and stress fields. However, the method 
requires a finite element mesh that is symmetric about the crack plane. 
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APPENDIX: NUMERICAL IMPLEMENTATION OF THE EDI METHOD 


This appendix presents the numerical implementation of the EDI method for 
a finite element analysis with isoparametric elements. The procedure presented in 
this appendix is similar to that of the 2-D analysis in reference 23 and is applicable 
to both 8-noded and 20-noded 3-D isoparametric elements. For the purpose of 
illustration the 20-noded element is used. 

A typical finite element model around the crack front is shown in Figure 
17(a). The shaded region represents a typical domain surrounding the crack front. 
Although, no restriction was imposed on the number of elements in the domain 
either in Xj — or X2 — directions, one ring of elements in xi — and X 2 — directions were 
used to explain the procedure. The procedure for computing J Xk is presented but 
can be easily extended to the 7/jj computation as well. 

The total 7-integral ( J Xh ) is equivalent to sum of the domain integrals con- 
tributed by the elements in the shaded region shown in Figure 17. 

( J*k)domain = 

«=1 

where J Xh is the volume integral over the i th element in the shaded region and 
N e is the number of elements enclosed in the domain. 

For isoparametric finite elements, the displacements within the element are 
defined by the shape functions Nj and the nodal displacements (u a )j. 

u a = (^2) 

where Nj = Nj(£, t), £) and £, 77, £ axe the coordinates of the parent element. 
The index j takes the value 1 to N n (N n is the number of nodes per element; N n 
= 8 for eight-node linear element and N n = 20 for twenty-node quadratic element) 
and a takes the value 1, 2, and 3 corresponding to the displacements in x 1 — , x 2 — , 
and X3— directions, respectively. 

The volume integral J Xh . of the i th element (Eq. (20)) is computed using 
Gaussian quadratures as 
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where Mq is the number of Gaussian quadrature points used in each direction £, rj, 
and £ and w m , w n , and w p are the Gaussian weights and del [J] is the determinant 
of the Jacobian matrix [J] defined by 


[J] = 


BN, BN, 0 %. 

~w ~st ' ’ ‘ ~~h 


BN, BN, 

8-n dn 


Bn 
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(A 4) 


Most of the quantities necessary for Equations (20-23), W , W 111 , {<r}, { 03 }i 
and [<t], are readily known in terms of the nodal displacements of the element. 
But, computation of the terms 5, {S’ 1 }, and needs special attention and is 
discussed below. Evaluation of the derivative matrices{cJ Cfc } and {u' 3 '} is same as 
that for hence, they are not discussed. 


5-Functions 

As mentioned previously, the 5-function is any arbitrary but continuous func- 
tion with a non-zero value (varying between 0 and 1) on the surface of the inner 
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tube (A c ) and a value of zero on the surface of the outer tube and at the two ends 
of the tube (see Fig. 1). The variation of 5-function over a typical i th element 
in the shaded region is shown in the Figure 17(b). The function is conveniently 
defined using the element shape functions as 


S(t, V, 0 = Nj Sj 

where j = 1 to N n and Sj is the nodal value of the 5-function at node j on the 
element. Different 5- functions can be defined by assigning 0, .5 or 1 to 5 j. For 
a typical element shown in the Figure 17(b), the 5-function is completely defined 
by specifying 52 = 5io = 1 and zero to all other nodes. This definition yields a 
5 -function having a parabolic variation along the crack front and a linear variation 
in the radial direction (Type I 5-function). 


Partial Derivatives of 5 

Once 5 is defined the partial derivatives of 5, 
computed using the isoparametric formulation as 


, GS . 


as ^ 

3TT 


W 

, es 

' dx 2 

> = [ J ]" 1 < 

as 

W 

&S 


as 

0*3 ' 


< &( j 


where [J] is the Jacobian matrix defined in Equation A4. 


377’ and 3^ can be 


(A6) 


Partial Derivatives of W 

The terms are computed by fitting a bilinear equation (in terms of the 
parent coordinates £, 77 , and £) to W , using the values at the 2 x 2 x 2 integration 
points and then taking derivatives with respect to x*. In reference 20, the integral 
was approximated by evaluating at the center of the element. A 
different approach is taken here. Because all the quantities are known at the 
integration points, the integration is carried out without further approximations 
of other terms in Equations (20) to (23). The values of the stresses are known to 
be more reliable at the 2x2x2 Gaussian points within the element (in comparison 
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to the nodal values). The stress- work density W is approximated in bilinear form 
as 

= ai + a 2 £ + a 3 77 + a 4 < + 05^77 

+ a 6 T)( + a 7 {( + a s (rjC (A7) 

Using the 2x2x2 Gaussian values of the stress-work density W, Equation (A7) 
is rewritten as 

W((,r),Q = [1 ( r, ( (, >,< « (VC HT]{W a } (AS) 

where 

111 11 1 11- 

-y/Z ~y/Z -V^ -v^ y/Z y/Z y/Z y/Z 
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3 3 _ 3 _3 _3 _3 3 3 ^ ’ 
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where Wi to Wvm are the values of W at the 2x2x2 Gaussian points shown 
in Figure 18. The partial derivatives 8X6 
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The derivatives 4^ can now be obtained as 
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where [J] is defined in Equation (A4). 

Similarly, derivatives of and can be obtained. All the necessary 

terms in Equations (20) to (24) are now known and, hence, the domain integrals 
for each element can be calculated. 

Computation of J-integral Along the Crack Front 

In a 3-D finite-element model of a cracked body, the crack front is divided into 
a number of To calculate the J-integral at each of the nodal points, for example 
at node i, consider the crack front segment between the nodes (t — 1) and (t + 1) 
(see Fig. 19). The 5-function will have a value of unity at node i and zero at 
nodes (i — 1) and ( i + 1). Since the 5-function is generated from the element shape 
functions, the 5-function is linear for 8-noded linear element and quadratic for the 
20-noded quadratic element (see Fig. 19). Utilizing the domain corresponding to 
the crack front segment ( i — 1) and ( i + 1) (see, for example Fig. 17 for a 20-node 
model), the J-integrals are calculated from Equations (15) and (16). The analysis 
is repeated at other nodal locations. If the first node is on the plane of symmetry 
of the model, the left-half (segment (t - 1) to i) of the 5-function is neglected. 
However, the accuracy of the J-integral is poor at nodes on the free surface of the 
model because of the well known boundary layer effect [38,39]. 
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Table 1. 5-functions and /-integrals. 


5-function type 

5-function 

/-integral 


Left segment 

Right segment 



C = x 3 /Aj 

II 

H 

CO 

> 

K> 


i 

(1-C 2 ) 

(1-C 2 ) 

2 A/3 

ii 

(1 + 0/2 

(1 - 0/2 

(Aj + A 2 V 2 

hi 

« 3 - 0/2 

(0 + 0/2 

(Ai + A2)/® 

IV 

« 2 - 0/2 

(0 + 0/2 

(Ai + A 2 )/6 

V 

(1+0/2 

(1-0/2 

(Ai + A2)/2 

VI 

(1+0/2 

(1-0/2 

(Ai + A2)/2 

8-noded element 

(1+0/2 

(1-0/2 

(Ai + A2)/2 


Table 2. Comparision of J Xl for various 5-functions and domains in a 
Kj = 1 stress field loaded specimen. 

5-function 

type 

D 2 

V'A. £ /(i - •'O 

Domains 

Z>3 D4 

D 5 

I 

1.0008 


1.0018 

II 

1.0008 


1.0018 

III 

1.0008 


1.0018 

IV 

1.0144 

1.0200 1.0085 

1.0062 

V 

1.0053 

1.0070 1.0037 

1.0032 

VI 

0.9870 

0.9870 1.0129 

0.9973 
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Table 3. Comparision of /-integral from EDI method with exact solu- 
tions for mixed-mode singular field loadings. 

(Type I 5-functions) 


K T 

Loading 

Kii 

Km 



Jx x 



Exact 

D 2 

EDI method 
Dz D4 

D 5 

1 

0 

0 

0.9100 

0.9107 

0.9105 

0.9112 

0.9116 

0 

1 

0 

0.9100 

0.9142 

0.9156 

0.9127 

0.9121 

0 

0 

1 

1.3000 

1.3015 

1.3020 

1.3011 

1.3010 

1 

1 

0 

1.8200 

1.7929 

1.7864 

1.8035 

1.8096 

0 

1 

1 

2.2100 

2.2158 

2.2176 

2.2138 

2.2131 

1 

0 

1 

2.2100 

2.2122 

2.2125 

2.2123 

2.2126 

1 

1 

1 

3.1200 

3.0945 

3.0885 

3.1046 

3.1106 


Table 4. Comparision of normalized J calculated from various domain 
definitions for M(T ) specimen. 

(W/a = 4, H/a = 8 , t/a = l,u = 0.3) 



E J Xl /{ 7 r <x 2 a 




EDI method 

VCCT method 

zjt Domain 


Domain type 



D a 

Db 

D c 

0.11 Dx 

1.558 

1.558 

1.558 1.575 

(interior layer) D 2 

1.551 

1.551 

1.571 

D 3 

1.544 

L559 

1.579 

D 4 

1.547 

1.572 

1.592 

D 5 

1.543 

1.582 

1.602 

0.495 I>i 

1.293 

1.293 

1.293 1.279 

(exterior layer) Dj 

1.380 

1.380 

1.296 

d 3 

1.468 

1.389 

1.305 

D 4 

1.498 

1.400 

1.316 

D s 

1.526 

1.411 

1.327 
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Figure 1.- Nomenclature and the domain description. 







Figure 2.- Symmetric and antisymmetric displacement components. 






finite-elements. 



XyU 
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Figure 6. Three-dimensional finite-element modeling of M(T) specimen. 



W/a 
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Figure 7. Variation of normalized J x along the crack front for a M(T) specimen. 
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Figure 8. Embedded penny-shaped crack in a circular rod subjected tension and torsion loadings. 
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Figure 9. Inclined penny-shaped crack embedded in a circular rod under tension 
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Figure 10. Finite-element model for an inclined embedded penny-shaped crack. 



1.00 
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Figure 1 1. Normalized J-distribution along the crack front of an embedded penny-shaped crack. 
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Figure 12(a). Normalized Jj, J n , and J m distribution along the crack front of an embedded inclined penny-shaped crack. 












Domain 
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Figure 12(b). Comparison of J m from direct and decomposition methods. 
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Figure 13. Normalized total J distribution along the crack front of an inclined semicircular surface crack. 
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Figure 14. 







a/c = 1/1.5 Green and Sneddon [35] 

t/a = 10 O EDI 

W/a = 10 
H/a = 222 


V 
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Figure 15. Application of EDI method to 8-node isoparametric element. 





c = 1 EDI 

t = 0.2 O VCCT [6] 

/a = 10 □ Newman & Raju [36] 
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Figure 16. Comparison of normalized J distribution for a semi-circular surface crack from EDI, VCCT, and force methods. 
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Figure 17.- Typical finite-element model and a domain description near the crack front. 
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Figure 18.- Nodes and 
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